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1. INTRODUCTION 

The Nigerian national power grid is powered from 29 generation stations having an installed 
capacity of 12,910.40MW but a relatively lower generation capacity of roughly 7,652.6MW. 19% of 
the installed capacity is provided by the more reliable and cheap hydropower plants. The three major 
hydropower plants in Nigeria are the Kainji hydroelectric power station (KHEPS), the Jebba hydroelectric 
power station (JHEPS) and the Shiroro hydroelectric power station [1-3]. 

The KHEPS located on 09°51'45"N, 04°36'48"E and JHEPS on 9°08'08"N, 4°47'16"E are in 
tandem and managed by the mainstream energy solution. They are built on the River Niger and operated in 
cascade. Unfortunately, the generation efficiency of JHEPS has been lower than expected, the daily average 
energy generated in 2017 was estimated at around 846 MW [4]. The two power stations operate in cascade 
but lack a control system regulating their operation. They are being managed by the experience and intuition 
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of the operators [2, 5]. From the operational report, there are occasions where some units at JHEPS are shut 
down if the release from KHEPS is low [6]. 

Based on the available information gathered, the problem of real-time optimal management of 
resources between KHEPS and JHEPS has not been solved. Despite all the research efforts made so far, 
operators still rely on intuitive water release rules to maximize power generation regularly. It 1s obvious that 
the operational procedure based on intuition and experience can leads to poor utilization of available energy 
resources [5, 7]. There are no proper scientifically motivated techniques for the management of resources 
between the two stations. Most models and optimization techniques earlier proposed are parameter 
optimization technique which is not sufficient for managing a dynamic system nonlinear system [8-11]. 

A difficult problem encountered by the plant operators occurs in the management of water flows 
between the two reservoirs. Experience from years of operations seems to support this approach given 
the fact that a necessary and enough condition for the proper operation of the turbo-alternators is that 
the operating net head must satisfy the requirements for acceptable energy conversion by the turbine. 
Therefore, the optimal control problem solved in this work is the determination of the best control vector and 
resulting state trajectories which minimizes certain performance index, subject to system constraints. In this 
case, the optimal control problem results, whereby a control signal is desired that will force the reservoir 
head at JHEPS to move from an initial point to the desired point in a finite time and subject to constraints 
imposed by the system dynamics. 

Unfortunately, many problems that are rooted in nonlinear optimal control theory do not have 
computable solutions or they have solutions that may be obtained only with a great deal of computing 
effort [12, 13]. The solution via analytical means also seems not feasible except by numerical means. Numerical 
solutions to optimal control problems are either through direct or indirect methods. In the direct methods, 
the infinitely dimensional state and controls are discretized. The indirect method applies calculus of variation to 
set up necessary conditions that must be satisfied by the optimal control. Calculus of variation, together with 
Pontryagin’s minimum principles are used to setup optimality conditions. These conditions produce optimal 
control canonical equations such that their solution ensures that an optimum point has been reached [14, 15]. 

The indirect approach leads to a nonlinear two-point boundary-value problem. The control task then 
reduces to the solution of a boundary value problem. There are different approaches with associated advantages 
and disadvantages. In all the solution techniques, an initial guess is used to obtain a solution in which one or 
more of the necessary optimality conditions are not satisfied. The solution is then used to adjust the initial guess 
to make the next solution come closer to satisfying all the necessary conditions. If the steps are repeated and 
the iterative procedure converges, the necessary conditions will eventually be reached [16-18]. 

Many authors had proposed methods of solving an optimal control problem, these methods can be in 
the form of nonlinear programming, shooting method, forward backward sweep, steepest descent, conjugate 
gradient, dynamic programming, the variation of externals, quasi-linearization, gradient projection, 
collocation, etc [19, 20]. There has been no perfect method as each has its own advantages and 
disadvantages. For example, the forward-backward sweep (FBS) works only if the Lipschitz constants for 
the state, costate and control variables are small enough or the time interval is very small. Likewise, 
the covergence of the shooting method depends on the numerical procedure and the initial data set, else there 
will be no solution [21]. 

The multiple shooting and parallel shooting techniques were earlier explored in [22-27] by resetting 
the problem and increasing the single initial value problem to a family of initial value problems configured so 
as to limit the effect of the growth of computational errors. The outcome resulted in a method that increased 
the number of guesses which were much fewer than the methods that depended on variational and other 
approximation methods that for the same accuracy may involve the solution of large linear or nonlinear 
equations that have dimensionality several orders of magnitude when compared with the corresponding 
shooting method. The method was applied successfully in the modeling of distributed parameter systems and 
proved to be very efficient, accurate and fast. The progressive domain expansion method (PDEM) proposed 
in this paper is another modification shooting method. It is less computational and feasible in solving 
the optimal control problem at JHEPS. 


2. RESEARCH METHOD 

The solution to an optimal control problem requires a proper model of the system dynamics in 
the form of a differential equation or difference equation. A suitable performance index with its associated 
constraints must be developed; there should also be a numerical technique for solving the model equation 
and a procedure for solving the resulting boundary value problem. The system can be described by 
the block diagram of Figure 1 and the dynamical model of (1). This model was earlier presented in [5]. 
The block diagram presents the reservoir dynamics and power generation as a function of the operating 
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head. As presented in [5], the electrical power form JHEPS is represented in (1) and (2), showing that it is 
a function of the operating head. 


P(t) = V2 npA, g /2h/2(t) (1) 
P(t) = K(h(t)) (2) 


where K represents a scalar valued function. 


VŽ npA, g !2 


© h h 
+ 
vo 


So a~ 
nj2g A ` Az 








JHEPS Q 
Upstream 








h (m): operating head 

Q (m3 /s): Inflow into the reservoir 

q (m/s): Inflow into the penstock 

0,(m? /s): Losses 

Q; (m? /s): Discharge through spill way 
A, (m*): Effective surface area 











Spillway 
function 







Q; 


A. (m*): Area of the scroll casing 
o(kg/m?*): Density of water 
g(ms*): Acceleration due to gravity 
n: Conversion efficiency 


JHEPS U, to Ua: The turbo-altemators 


discharge 


=q+ Qs 





Figure 1. Block diagram of the JHEPS reservoir dynamics 


Therefore, (3) presents the differential equation (model equation) that describes the operating 
head dynamics, 


MO = L nA, tA, [2gh(t) +A, OO) - OE — O) (3) 
h(t) = f (A(t), u(t)) (4) 


where f also represents a scalar valued function. 

The optimal control of the JHEPS devolves into the determination of best control vector 
u(t) € U(t) from KHEPS, which compels the dynamical system h(t) = f (h(t), u(t), t) at JHEPS to follow 
an optimal trajectories h*(t) that minimize specified performance indices. The system is nonlinear and exists 
in the continuous-time domain. The optimal control problem is the Lagrange problem. As earlier mentioned, 
while solving the most optimal control problems, the solution can be intractable by analytical methods; 
the engineers depend on numerical procedures. The numerical solution used in determining the optimal 
solution is the progressive expansion of domain method. Consequently, the problem addressed here adopts 
a performance index that includes penalties for deviation from a specified state and for deviations from some 
predefined value of the control variable. 


J = min fif (Kp (hŒ) — A(T)” + K, (ult) - u(r)’ ) at (5) 

J = min ff oh), uO), t) dt; (6) 
Subject to: 

h(t) = fh u(t), t) 

h(to) = ho (7) 
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h(t,;) = ACT) (8) 
j= tat 


where ¢ is a function, K, and K,, are non-zero positive constants. 
The (5) can be written as an augmented optimization problem by introducing a Lagrange multiplier 
A(t) as follows: 








J = S (PAOU, +20 (FROUD - hO)) 4 (9) 
applying the method of Pontryagin, a Hamiltonian function H (h(t), u(t), A(t), t) is defined as: 
H = p(h@), ult), t) +AM FW), ult), t) (10) 
jej (H = A(t)h(t)) dt (11) 
applying the Euler-Lagrange equation to the Hamiltonian gives: 
OF. Wis 
an(t) J. (t)| = 
A) =- (12) 
ðH ; 
aw 7 f AW), ue), t) — AE) = 0 (13) 
h(t) = f(A), u(t), t) (14) 
ðH 
dult) (a) 
therefore, 
H = K,(h(t) -ACT))? + Kalule) = UCT)? + AC) (rah + wuce)) 
A(t) = —njah@/2)(t) + p ult) (16) 
A(t) = —[2K,(h(t) — ACT)) — 2a (nao) (17) 
0 = 2K,(u(t) —u(T)) + pact) (18) 
and the boundary condition, 
A(to) = ho 
A(t,) =1 (19) 


From the result above, the determination of the optimal control requires the solution of a two-point 
boundary value problem (TPBVP) since the initial conditions of the system are specified at the initial time 
and the value of the Lagrange multiplier or co-state is specified the terminal point. Meanwhile, the unknown 
control is related to the state and co-state through the optimality condition. From (18), w(t) can be expressed as: 


u(t) = -H + u(r) (20) 


substituting (20) into (16) leads to two sets of first-order differential equations with split boundary conditions. 
h(t) = -njah e) — WES + u(r) (21) 
A(t) = -2K (ACE) — A(T) +24 (nande) 
h(tọ) = họ and h(t;) = unknown 


A(to) = unknown and Afte) =i (22) 
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2.1. Computation of the optimal control by a progressive domain expansion method (PDEM) 

In (21) and (22) provide the basis for the computation of optimal control that consists of two 
nonlinear differential equations with split boundary conditions. The state equation is specified at the initial 
time while the co-state equation has its condition specified at the final time. The two-point boundary value 
problem (TPBVP) has attracted considerable attention in the past century. While many methods are found 
in the literature, one method that was often proposed but seldom used because of computational difficulties is 
the shooting technique. It is attractive because it involves guessing values for the missing initial conditions 
which in principle ought to be determined by one of several possibilities, the nature of the co-state equations 
lead to a rapid growth of the initial value problem that the computed values soon lose relationship with 
the problem since errors in computation are exponentially amplified by the system. 

The PDEM is another modification proposed and pre-tested by [26], whereby instead of partitioning 
the domain [0, T], the final domain boundary is adjusted in such a manner that the initial guess results still 
retain a semblance with the original problem and the growing equations are bound so that the correct initial 
value problem is solved assuming the pseudo-domain [0,7;] where Tg is determined on the fly. In the next 
iteration, the initial guess for the missing boundary condition assumes the value that would have solved 
the problem for the pseudo-domain. Meanwhile, the problem is solved beyond Tę again until the growth 
begins to cause concern. The new 7;, 1s thus defined and a new correction made so that the corresponding 
problem is solved. This process is repeated until Tę = Tp in which case the value of the initial guess 
converges to the correct initial guess for the problem. The procedure is presented in the flowchart of Figure 2. 


fs (h Ct), A® Ct), t) 





fi(hP O, AO), o) _ | —n ah (t)'/2 — 2K, +yuu(T) 


pwea™(t) | 


1 : 
—2Khn (hC) = h(T)) + = nja Agh (t) 2 
Let k=0 


h(t) = ho 
Guess value for A‘*=°) (to) 


Numerically integrate the state and 
co-state equation from ty > t; < ty 


Yes 





No 


Perturb 4%?) (t6) 










If 
maxfh (t), A (t7)} <M 
M is within computer 
precision. 


Obtain the new state and 
co-state h™ (tj) and A™(t;) 











A’ Cto) = AM Cto) 
Obtain the optimal control 
u(t) and the optimal state 
trajectories h(t) 
uË Ct) 
2K,, 


fort = to>ty; 





If 
Ja(t,) — 1] < 10™” 
nis a positive integer 


Yes 


TON i . No 
Perturb A‘ (to) with 5A(ty) 
and obtain dA(t,) 


ak+D(t,,) = a(t) z (A (t,) — 1) (Sin ) 


u(t)= — 





+u(T); 


=} 





End 


Figure 2. Flow chart for the solution of optimal control canonical equations by a PDEM 
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3. RESULTS AND ANALYSIS 

It is better to define a presentation format that reflects important features of the trajectory, many of 
which are deduced but not definitive by themselves in assessing a given trajectory. This leads to a very 
important matter of how to use optimal control for operational purposes. The notations for specifying 
operating conditions were formulated as follows: (number of operating machines, starting head (m), duration 
(hr), constraints on maximum inflow). This notation would be used in the presentation of results. 


3.1. Case 1: (5,25.8,24,u(T) unconstrained) 

Figure 3 presents the result for an indirect optimal control with the operational conditions of 
the machines and the initial value of the head at JHEPS specified as (5,25.8,24,u(T) unconstrained). 
The operating condition implies that given that 5 turbo-alternators are running, while the operating head is 
25.8 m. It is desired that the operating head increases to the nominal value 26.1m in 24 hr (86400s). 
The control problem is the determination of the inflow required to achieve this stated objective. The optimal 
control system generated a control law of (23); 


u(t) = —1x107*'t? + 2x 107°t* — 0.1182t + 4715.9 (23) 


The control started from a value around 4605 m?/s at tọ and decreases gradually to zero at tr. Because 
the maximum control is not constrained, the trajectory of the operating head could not reach the terminal 
value h(T) = 26.1, as h(t; ) = 25.94m. The operating head also rose to a peak value and decreases later. 
Hence this result is not satisfactory, the algorithm has to be modified. 


3.2. Case 2: (3,25.8,24,u(T) unconstrained) 

Case 2 considered a situation where the number of units in operation reduces to 3 machines and 
the operating head 25.8 m. A similar result to that of Figure 3 was observed and presented in Figure 4. It can 
be concluded that to use PDEM, a constraint in the form of penalty on the terminal control u(tr) can be 
incorporated into the algorithm. 


u(t) = -9x10°77t? + 2x 10°°t? — 0.1165t + 4143.5 (24) 
== Inflow (m3/s) —@— Head(m) —®— Inflow (m3/s) —@®— Head(m) 
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3.3. Case 3: (5,25.8,1,u(T) = 1800 m4) 

Case 3 presents a similar condition to that of case 1 but with a constraint on the final control. This 
implies that the final control cannot decrease to zero, but a finite value specified as u(T). The procedure for 
the selection of suitable value for w(T) can be found in [5]. Figure 5 shows a better result of which the head 
moves from an initial value h(0) = 25.8 m, to a final value h(T) = 26.1 m without an overshoot, as 
experienced in the direct optimal control. The optimal control defined as (25) starts around 4327.4 m? /s and 
ends at 1800 m?/s. 


u(t) = 5x 107°7t? — 0.0702t + 4327.4 (25) 
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3.4. Case 4: (5,25.8,2,u(T) = 1800 m°*%) 

A situation can be considered where the time limit for which optimal control is required increases to 
two days (48 hr.). It can be observed from Figure 6 that the PDEM could determine the solution as well 
except. In case 3, u(0) started with 4327.4 m? /s but reduces slightly to 4270.1 m?/s in case 4, this implies 
that the time does not have significant effect on the starting control u(0). The control law, in this case, is 
expressed as (26). 


u(t) = —2 x 10712t3 + 6x 107°7t? — 0.0689t + 4270.1 (26) 
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3.5. Approximate optimal control: (5, 25.8, 1,u(T) = 1800 m?/s) 

Consider the realization of the physical controller or the release of inflow from KHEPS, the infinite 
dimensional optimal control generated in case 1 to case 4 may not be easy to implement. An approximate 
optimal control may be necessary such that the inflow is released systematically withing every 6 hrs. 
This was considered for case 1, such that the average control with every 6 hrs was determined, the control 
is presented in Figure 7, where ul = 3648.20 m?/s,u2 = 2544.90 m?/s,u3 = 2040.46m?/s, 
u4 = 1837.07 m?/s. To show the performance of this approximation, the operating head trajectory 
generated is presented in Figure 8 to be satisfactory. 
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4. CONCLUSION 

This paper proposed a progressive expansion of domain technique as a numerical approximation to 
the solution of an optimal control problem involving the regulation of the operating head of JHEPS. 
The solution was necessary knowing that the resulting two-point boundary value problem imposes 
a limitation in the use of indirect optimal control. From the result, the technique does not require 
a sophisticated initial guess like the normal shooting techniques and it converges faster than results expected 
from the nonlinear programming approach. The algorithm is therefore recommended for use in system 
studies, management and complete design of an electronic controller for the station. 
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